
# step 12 Finding differentially expressed features (cluster biomarkers)
#1 find all markers of cluster 3
cluster3.markers <- FindMarkers(pbmc, ident.1 = 3, min.pct = 0.25)
#  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
head(cluster3.markers, n = 5)


#2
cluster4_02.markers <- FindMarkers(pbmc, ident.1 = 4, ident.2 = c(0, 2), min.pct = 0.25)
#  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
head(cluster4_02.markers, n = 10)




$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -n "^\s*FindMarkers" 2>/dev/null --color=auto
./seurat-4.1.0/R/differential_expression.R:113:        FindMarkers(
./seurat-4.1.0/R/differential_expression.R:492:FindMarkers.default <- function(
./seurat-4.1.0/R/differential_expression.R:611:FindMarkers.Assay <- function(
./seurat-4.1.0/R/differential_expression.R:690:FindMarkers.SCTAssay <- function(
./seurat-4.1.0/R/differential_expression.R:793:FindMarkers.DimReduc <- function(
./seurat-4.1.0/R/differential_expression.R:909:FindMarkers.Seurat <- function(
./seurat-4.1.0/R/generics.R:143:FindMarkers <- function(object, ...) {





() 泛型函数
../seurat-4.1.0/R/generics.R:143:FindMarkers <- function(object, ...) {

#' Gene expression markers of identity classes
#'
#' Finds markers (differentially expressed genes) for identity classes
#'
#' @param object An object
#' @param ... Arguments passed to other methods and to specific DE methods

#' @return data.frame with a ranked list of putative markers as rows, and associated
#' statistics as columns (p-values, ROC score, etc., depending on the test used (\code{test.use})). The following columns are always present:
#' \itemize{
#'   \item \code{avg_logFC}: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group
这里v4版本已经是 avg_log2FC： 2组间平均表达量的 log变化倍数。
#'   \item \code{pct.1}: The percentage of cells where the gene is detected in the first group 第一组中细胞的表达的百分比
#'   \item \code{pct.2}: The percentage of cells where the gene is detected in the second group 第一组中细胞的表达的百分比
#'   \item \code{p_val_adj}: Adjusted p-value, based on bonferroni correction using all genes in the dataset 使用 bonferroni 方法矫正后的p值
#' }
#'
#' @details p-value adjustment is performed using bonferroni correction based on
#' the total number of genes in the dataset. Other correction methods are not
#' recommended, as Seurat pre-filters genes using the arguments above, reducing
#' the number of tests performed. Lastly, as Aaron Lun has pointed out, p-values
#' should be interpreted cautiously, as the genes used for clustering are the
#' same genes tested for differential expression.
p值矫正使用的是严格的 bonferroni，不建议使用其他矫正方法，因为事先过滤基因了，导致检验数少了很多。
最后，像 Aaron Lun 指出的，p值要谨慎解释，因为分组和差异表达使用的同样的基因。

#' 
#' @references McDavid A, Finak G, Chattopadyay PK, et al. Data exploration,
#' quality control and testing in single-cell qPCR-based gene expression experiments.
#' Bioinformatics. 2013;29(4):461-467. doi:10.1093/bioinformatics/bts714 
检验

#' @references Trapnell C, et al. The dynamics and regulators of cell fate
#' decisions are revealed by pseudotemporal ordering of single cells. Nature
#' Biotechnology volume 32, pages 381-386 (2014)
pseudotemporal

#' @references Andrew McDavid, Greg Finak and Masanao Yajima (2017). MAST: Model-based
#' Analysis of Single Cell Transcriptomics. R package version 1.2.1.
#' https://github.com/RGLab/MAST/
另一个包 MAST

#' @references Love MI, Huber W and Anders S (2014). "Moderated estimation of
#' fold change and dispersion for RNA-seq data with DESeq2." Genome Biology.
#' https://bioconductor.org/packages/release/bioc/html/DESeq2.html
DESeq2 分析RNA-seq的倍数变化和分布

#'
#' @export
#'
#' @examples
#' data("pbmc_small")
#' # Find markers for cluster 2
#' markers <- FindMarkers(object = pbmc_small, ident.1 = 2)
#' head(x = markers)
#'
#' # Take all cells in cluster 2, and find markers that separate cells in the 'g1' group (metadata
#' # variable 'group')
#' markers <- FindMarkers(pbmc_small, ident.1 = "g1", group.by = 'groups', subset.ident = "2")
#' head(x = markers)
#'
#' # Pass 'clustertree' or an object of class phylo to ident.1 and
#' # a node to ident.2 as a replacement for FindMarkersNode
#' if (requireNamespace("ape", quietly = TRUE)) {
#'   pbmc_small <- BuildClusterTree(object = pbmc_small)
#'   markers <- FindMarkers(object = pbmc_small, ident.1 = 'clustertree', ident.2 = 5)
#'   head(x = markers)
#' }
#'
#' @rdname FindMarkers
#' @export FindMarkers
#'
#' @aliases FindMarkersNode
#' @seealso \code{FoldChange}
#'
FindMarkers <- function(object, ...) {
  UseMethod(generic = 'FindMarkers', object = object)
}







()
../seurat-4.1.0/R/differential_expression.R:909:FindMarkers.Seurat <- function(

#' @param ident.1 Identity class to define markers for; pass an object of class
#' \code{phylo} or 'clustertree' to find markers for a node in a cluster tree;
#' passing 'clustertree' requires \code{\link{BuildClusterTree}} to have been run
参数 ident.1 可以是要查找 marker 的类。 或者是 phylo 或 clustertree 类聚类树的一个节点。


#' @param ident.2 A second identity class for comparison; 
#if \code{NULL}, use all other cells for comparison; 
#if an object of class \code{phylo} or 'clustertree' is passed to \code{ident.1}, must pass a node to find markers for
参数 ident.2，另一个要比较的身份类。
- 如果为空，表示使用其他细胞。
- 如果 ident.1 传入的是一个 phylo 或 clstertree 类型，则ident.2必须传入一个要找marker的节点。 


#' @param reduction Reduction to use in differential expression testing - will test for DE on cell embeddings
reduction: 差异表达检验使用的降维方法 - 在细胞嵌入上做DE检验


#' @param group.by Regroup cells into a different identity class prior to performing differential expression (see example)
差异表达分析前，对细胞分组，按照 group.by


#' @param subset.ident Subset a particular identity class prior to regrouping. Only relevant if group.by is set (see example)
分组前，取子集。只有设置 group.by 时有效。


#' @param assay Assay to use in differential expression testing
使用哪个 assay做差异分析检验


#' @param slot Slot to pull data from; 
slot: 从哪个slot获取数据
# note that if \code{test.use} is "negbinom", "poisson", or "DESeq2", \code{slot} will be set to "counts"
注意: 如果检验方法是 "negbinom", "poisson", or "DESeq2"，则设置从 counts 获取矩阵。



#' @param mean.fxn Function to use for fold change or average difference calculation.
#' If NULL, the appropriate function will be chose according to the slot used
mean.fxn: 用于计算 倍数变化 或 平均差异的函数。如果为 NULL，会根据slot选择合适的函数。



#' @param fc.name Name of the fold change, average difference, or custom function column in the output data.frame. 
#If NULL, the fold change column will be named according to the logarithm base (eg, "avg_log2FC"), 
# or if using the scale.data slot "avg_diff".
fc.name: 输出数据框中 倍数变化、平均差异、或自定义函数列的名字。
如果为NULL，倍数变化列名字将会根据log底命名，比如 avg_log2FC；如果使用 scale.data slot 则命名为 avg_diff


#' @param base The base with respect to which logarithms are computed.
base: 对数底用哪个？ 默认是2

#'
#' @rdname FindMarkers
#' @concept differential_expression
#' @export
#' @method FindMarkers Seurat
#'
FindMarkers.Seurat <- function(
  object,
  ident.1 = NULL,
  ident.2 = NULL,
  group.by = NULL,
  subset.ident = NULL,
  assay = NULL,
  slot = 'data',
  reduction = NULL,
  features = NULL,
  logfc.threshold = 0.25,
  test.use = "wilcox",    # 默认使用 wilcox 检验
  min.pct = 0.1,
  min.diff.pct = -Inf,
  verbose = TRUE,
  only.pos = FALSE,
  max.cells.per.ident = Inf,
  random.seed = 1,
  latent.vars = NULL,
  min.cells.feature = 3,
  min.cells.group = 3,
  pseudocount.use = 1,
  mean.fxn = NULL,
  fc.name = NULL,
  base = 2,
  densify = FALSE,
  ...
) {
  #(A1) 如果 group.by 非空
  if (!is.null(x = group.by)) {
  	# 如果 subset.ident 非空
    if (!is.null(x = subset.ident)) {
      # 先取子集
      object <- subset(x = object, idents = subset.ident)
    }
    # 然后改为使用 group.by 作为 Idents
    Idents(object = object) <- group.by # getter见 解析9-2.2, setter 见本文
  }
  # 关键步骤。小结: 先取子集，然后改为 group.by 设为 slot active.ident


  #(A2) 如果同时指定 assay 和 reduction，则报错: 只能指定一个
  if (!is.null(x = assay) && !is.null(x = reduction)) {
    stop("Please only specify either assay or reduction.")
  }

  #(A3) 如果 ident.1 长度为0，则报错
  if (length(x = ident.1) == 0) {
    stop("At least 1 ident must be specified in `ident.1`")
  }

  #(A4) 选择使用的数据 # select which data to use
  # 如果 reduction 为空，则使用 assay
  if (is.null(x = reduction)) {
  	# > DefaultAssay(pbmc)
    # [1] "RNA"
    assay <- assay %||% DefaultAssay(object = object) #assay 默认值是 "RNA"
    data.use <- object[[assay]] #获得的是 assay 对象
    # > class( pbmc[[DefaultAssay(pbmc)]] )
    #[1] "Assay"
    cellnames.use <-  colnames(x = data.use) #获得列名: cell id
  } else {
  	# 如果 reduction 非空，则使用 pca、umap 数据?
  	# > names(pbmc@reductions)
    #[1] "pca"  "umap" "tsne"
    data.use <- object[[reduction]]
    cellnames.use <- rownames(x = data.use)
  }

  #(A5) 新函数: 立刻追踪该函数
  # 根据分类获得 cell id，返回一个list
  cells <- IdentsToCells(
    object = object,
    ident.1 = ident.1,
    ident.2 = ident.2,
    cellnames.use = cellnames.use
  )

  #(A6) 如果潜在变量非空，则获取
  # fetch latent.vars
  if (!is.null(x = latent.vars)) {
    latent.vars <- FetchData( #该函数见 解析 9-2.1
      object = object,
      vars = latent.vars,
      cells = c(cells$cells.1, cells$cells.2) #限定这几个细胞的属性
    )
  }

  # (A7) 检查标准化方法
  # check normalization method
  norm.command <- paste0("NormalizeData.", assay) #默认是 NormalizeData.RNA

  # (A8) 如果该方法在历史命令中，且 reduction 为空
  # > Command(pbmc)
  #[1] "NormalizeData.RNA"        "FindVariableFeatures.RNA" "ScaleData.RNA"            "RunPCA.RNA"              
  #[5] "FindNeighbors.RNA.pca"    "FindClusters"             "RunUMAP.RNA.pca"          "RunTSNE" 
  if (norm.command %in% Command(object = object) && is.null(x = reduction)) {
  	#(B1) 则标准化方法为: LogNormalize
  	#> Command(pbmc, "NormalizeData.RNA")
	#Command: NormalizeData(sce)
	#Time: 2022-01-09 16:25:11
	#assay : RNA 
	#normalization.method : LogNormalize 
	#scale.factor : 10000 
	#margin : 1 
	#verbose : TRUE
    norm.method <- Command(
      object = object,
      command = norm.command,
      value = "normalization.method"
    )
    #(B2) 如果标准化方法不是 LogNormalize
    if (norm.method != "LogNormalize") {
      # 则 求平均值的方法是: log2(行均值+1)
      mean.fxn <- function(x) {
        return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
      }
    }
  }
  # 如果标准化方法是 LogNormalize，则 mean.fxn 为空


  # (A9) 结果推给 .assay
  de.results <- FindMarkers(
    object = data.use,
    slot = slot,
    cells.1 = cells$cells.1, #输出的是 cell id list
    cells.2 = cells$cells.2,
    features = features,
    logfc.threshold = logfc.threshold,
    test.use = test.use,
    min.pct = min.pct,
    min.diff.pct = min.diff.pct,
    verbose = verbose,
    only.pos = only.pos,
    max.cells.per.ident = max.cells.per.ident,
    random.seed = random.seed,
    latent.vars = latent.vars,
    min.cells.feature = min.cells.feature,
    min.cells.group = min.cells.group,
    pseudocount.use = pseudocount.use,
    mean.fxn = mean.fxn,
    base = base,
    fc.name = fc.name,
    densify = densify,
    ...
  )


  return(de.results)
}






() IdentsToCells
$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -in "IdentsToCells" 2>/dev/null --color=auto
./seurat-4.1.0/R/differential_expression.R:1603:IdentsToCells <- function(

FindMarkers.Seurat 和 FoldChange.Seurat 的辅助函数.

把 idents 转为 细胞，就是返回一个list，名字是 idents，value是属于该idents的cell id.


# Helper function for FindMarkers.Seurat and FoldChange.Seurat
# Convert idents to cells
#
#' @importFrom methods is
#
IdentsToCells <- function(
  object,
  ident.1,
  ident.2,
  cellnames.use
) {
  #(A1) 如果 ident.1 为空，则报错
  if (is.null(x = ident.1)) {
    stop("Please provide ident.1")
    #(B1) 如果 ident.1 长度是1，且第一个值是字符串 "clustertree"，或者 ident.1 是类 'phylo' 的实例
  } else if ((length(x = ident.1) == 1 && ident.1[1] == 'clustertree') || is(object = ident.1, class2 = 'phylo')) {
  	# (C1) 如果 ident.2 为空，报错
    if (is.null(x = ident.2)) {
      stop("Please pass a node to 'ident.2' to run FindMarkers on a tree") #对树进行搜索时，1为树本身，2为要找的cluster
    }

    # (C2)如果 ident.1 是 系统发育树，则 tree = ident.1
    tree <- if (is(object = ident.1, class2 = 'phylo')) {
      ident.1
    } else {
    # 如果 ident.1 不是 系统发育树，则实用工具构建树 // 这部分跳过，水有点深，对树也不了解
      Tool(object = object, slot = 'BuildClusterTree')
    }

    # (C3)如果树为空，则报错
    if (is.null(x = tree)) {
      stop("Please run 'BuildClusterTree' or pass an object of class 'phylo' as 'ident.1'")
    }

    # (C4) 分别取出左枝和右枝
    ident.1 <- tree$tip.label[GetLeftDescendants(tree = tree, node = ident.2)]
    ident.2 <- tree$tip.label[GetRightDescendants(tree = tree, node = ident.2)]
  }


  #(A2) 如果 ident.1 长度大于1，且至少一个 ident.1 在 cellnames.use 中
  # 看来可以传入 cell id list 作为 ident.1
  if (length(x = as.vector(x = ident.1)) > 1 &&
      any(as.character(x = ident.1) %in% cellnames.use)) {
  	# 则找到 ident.1 中不在 cellnames.use 中的cell id
    bad.cells <- cellnames.use[which(x = !as.character(x = ident.1) %in% cellnames.use)]
    if (length(x = bad.cells) > 0) {
      # 如果有不在的，则报错
      stop(paste0("The following cell names provided to ident.1 are not present in the object: ", paste(bad.cells, collapse = ", ")))
    }
  } else {
  	# 否则，直接获取属于 ident.1 的cell id
    ident.1 <- WhichCells(object = object, idents = ident.1)
  }


  #(A3) 处理 ident.2: 如果为空，则使用其他全部细胞
  # if NULL for ident.2, use all other cells
  # (B1)如果长度大于1，且至少有一个细胞在 cellnames.use 中
  if (length(x = as.vector(x = ident.2)) > 1 &&
      any(as.character(x = ident.2) %in% cellnames.use)) {
  	# 则找 ident.2 不在 cellnames.use 的细胞，如果有，则报错
    bad.cells <- cellnames.use[which(!as.character(x = ident.2) %in% cellnames.use)]
    if (length(x = bad.cells) > 0) {
      stop(paste0("The following cell names provided to ident.2 are not present in the object: ", paste(bad.cells, collapse = ", ")))
    }
  # (B2) 否则，就是2个情况: ==1 和 ==0
  } else {
  	# 如果为空，就是==0，则使用其余全部细胞
    if (is.null(x = ident.2)) {
      ident.2 <- setdiff(x = cellnames.use, y = ident.1)
    } else {
    # 否则，就剩下 ==1，则取出该类的细胞
      ident.2 <- WhichCells(object = object, idents = ident.2)
    }
  }

  #(A4) 返回细胞list: key为分类，value为该分类下的 cell id
  return(list(cells.1 = ident.1, cells.2 = ident.2))
}








()FindMarkers.Assay()
../seurat-4.1.0/R/differential_expression.R:611:FindMarkers.Assay <- function(

#' @rdname FindMarkers
#' @concept differential_expression
#' @export
#' @method FindMarkers Assay
#'
FindMarkers.Assay <- function(
  object,
  slot = "data",
  cells.1 = NULL,
  cells.2 = NULL,
  features = NULL,
  logfc.threshold = 0.25,
  test.use = 'wilcox',
  min.pct = 0.1,
  min.diff.pct = -Inf,
  verbose = TRUE,
  only.pos = FALSE,
  max.cells.per.ident = Inf,
  random.seed = 1,
  latent.vars = NULL,
  min.cells.feature = 3,
  min.cells.group = 3,
  pseudocount.use = 1,
  mean.fxn = NULL,
  fc.name = NULL,
  base = 2,
  densify = FALSE,
  ...
) {
  #(A1) 根据统计检验的方法，决定使用 counts 还是标准化后的 data
  # 我们使用默认值 wilcox，则这里 data.slot="data"
  data.slot <- ifelse(
    test = test.use %in% DEmethods_counts(),
    yes = 'counts',
    no = slot
  )
# > Seurat:::DEmethods_counts
# function () 
# {
#     c("negbinom", "poisson", "DESeq2") #这3个需要原始counts
# }
# <bytecode: 0x55b94f4a0738>
# <environment: namespace:Seurat>
# > Seurat:::DEmethods_counts()
# [1] "negbinom" "poisson"  "DESeq2"  


  #(A2) 相当于直接使用 obj@data
  data.use <-  GetAssayData(object = object, slot = data.slot)
 
# > obj=pbmc@assays$RNA
# > class(obj)
# [1] "Assay"
# attr(,"package")
# [1] "SeuratObject"

# > data.use=obj@data #函数实际执行的句子

# > data.use2 = GetAssayData(object = obj, slot = "data")
#> unlist(data.use - data.use2) |> as.numeric() |> table()
#       0 
#36177532 

  #(A3) 默认 data.slot =="data"，所以 counts=numeric()
  counts <- switch(
    EXPR = data.slot,

    # 如果想要 scale.data，这里获取 counts?
    'scale.data' = GetAssayData(object = object, slot = "counts"),
    numeric()
  )


  #(A4) 新函数: 立刻查看。返回3列df: log2FC, pct.1, pct.2 
  fc.results <- FoldChange(
    object = object,
    slot = data.slot,
    cells.1 = cells.1,
    cells.2 = cells.2,
    features = features,
    pseudocount.use = pseudocount.use,
    mean.fxn = mean.fxn,
    fc.name = fc.name,
    base = base
  )

  #(A5) 计算任务推给 FindMarkers.default: p值还没有计算
  de.results <- FindMarkers(
    object = data.use,
    slot = data.slot,
    counts = counts,
    cells.1 = cells.1,
    cells.2 = cells.2,
    features = features,
    logfc.threshold = logfc.threshold,
    test.use = test.use,
    min.pct = min.pct,
    min.diff.pct = min.diff.pct,
    verbose = verbose,
    only.pos = only.pos,
    max.cells.per.ident = max.cells.per.ident,
    random.seed = random.seed,
    latent.vars = latent.vars,
    min.cells.feature = min.cells.feature,
    min.cells.group = min.cells.group,
    pseudocount.use = pseudocount.use,
    fc.results = fc.results,
    densify = densify,
    ...
  )


  return(de.results)
}














()FoldChange 泛型函数
FindMarkers.Assay() 中调用的应该是 FoldChange.Assay()。
不过，我们把同名函数都看一眼吧。

$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -n "^\s*FoldChange" 2>/dev/null --color=auto

./seurat-4.1.0/R/differential_expression.R:1093:  FoldChange( #调用 defautl 函数
#如果是 pca 或者 umap数据。 //todo 略过
./seurat-4.1.0/R/differential_expression.R:1108:FoldChange.DimReduc <- function( 


泛型函数:
./seurat-4.1.0/R/generics.R:242:FoldChange <- function(object, ...) {

计算不同类别之间的: 每个基因的 log 倍数变化，和细胞表达百分比。

如果使用 scale.data 或者 降维数据，返回平均差异，列名为 avg_diff;
否则返回 log 倍数变化，列名为 avg_log2_FC

实例
# > FoldChange(pbmc, ident.1 = 1) |> head(n=3)
#                avg_log2FC pct.1 pct.2
# AL627309.1     0.05977664 0.012 0.001
# AP006222.2     0.01078354 0.002 0.001
# RP11-206L10.2 -0.01000233 0.000 0.002

#' Fold Change
#'
#' Calculate log fold change and percentage of cells expressing each feature
#' for different identity classes.
#'
#' If the slot is \code{scale.data} or a reduction is specified, average difference
#' is returned instead of log fold change and the column is named "avg_diff".
#' Otherwise, log2 fold change is returned with column named "avg_log2_FC".
#'
#' @examples
#' data("pbmc_small")
#' FoldChange(pbmc_small, ident.1 = 1)
#'
#' @param object A Seurat object
#' @param ... Arguments passed to other methods
#' @rdname FoldChange
#' @export FoldChange
#' @return Returns a data.frame
#' @seealso \code{FindMarkers}
FoldChange <- function(object, ...) {
  UseMethod(generic = 'FoldChange', object = object)
}






() FoldChange.Seurat()
./seurat-4.1.0/R/differential_expression.R:1156:FoldChange.Seurat <- function(



#' @param ident.1 Identity class to calculate fold change for; 
# pass an object of class \code{phylo} or 'clustertree' to calculate fold change 
# for a node in a cluster tree;
#' passing 'clustertree' requires \code{\link{BuildClusterTree}} to have been run
ident.1: 计算倍数变化的类别名称1
- 对于树中的节点，输入一个 phylo 对象，或者 字符串 'clustertree'，后者需要运行 BuildClusterTree

#' @param ident.2 A second identity class for comparison; if \code{NULL},
#' use all other cells for comparison; if an object of class \code{phylo} or
#' 'clustertree' is passed to \code{ident.1}, must pass a node to calculate fold change for
ident.2: 比较的第二个参数；
如果为空，则使用其余所有细胞；
如果第一个参数是 phylo 对象 或者 'clustertree'，则第二个参数必须是一个节点。


#' @param reduction Reduction to use - will calculate average difference on cell embeddings
reduction: 计算细胞嵌入点的平均差异。

#' @param group.by Regroup cells into a different identity class prior to
#' calculating fold change (see example in \code{\link{FindMarkers}})
group.by: 计算倍数变化前，重新把细胞分类

#' @param subset.ident Subset a particular identity class prior to regrouping.
#' Only relevant if group.by is set (see example in \code{\link{FindMarkers}})
subset.ident: 使用 group.by 前，取子集。
group.by 参数设置时，该参数才会有意义。

#' @param assay Assay to use in fold change calculation | 计算倍数变化的 Assay 类
#' @param slot Slot to pull data from | 从哪个slot获取数据
#' @param pseudocount.use Pseudocount to add to averaged expression values when
#' calculating logFC. 1 by default.
pseudocount.use: 计算log倍数变化时，防止log无意义，对counts默认加1

#' @param mean.fxn Function to use for fold change or average difference calculation
mean.fxn: 计算倍数变化或 平均差异的函数

#' @param base The base with respect to which logarithms are computed.
base: 倍数的底，默认为2

#' @param fc.name Name of the fold change, average difference, or custom function column
#' in the output data.frame
fc.name: 倍数变化的列名。

#'
#' @rdname FoldChange
#' @concept differential_expression
#' @export
#' @method FoldChange Seurat
FoldChange.Seurat <- function(
  object,
  ident.1 = NULL,
  ident.2 = NULL,
  group.by = NULL,
  subset.ident = NULL,
  assay = NULL,
  slot = 'data',
  reduction = NULL,
  features = NULL,
  pseudocount.use = 1,
  mean.fxn = NULL,
  base = 2,
  fc.name = NULL,
  ...
) {
  #(A1) 如果 group.by 非空
  if (!is.null(x = group.by)) {
  	# 如果 subset.ident 非空，则取子集
    if (!is.null(x = subset.ident)) {
      object <- subset(x = object, idents = subset.ident)
    }
    # 把 group.by 设置为 新分类: active.ident slot
    Idents(object = object) <- group.by
  }

  #(A2) 如果assay和 reduction 同时非空，则报错: 只能设置一个
  if (!is.null(x = assay) && !is.null(x = reduction)) {
    stop("Please only specify either assay or reduction.")
  }

  #(A3) 使用哪个数据? select which data to use
  # 如果 reduction 为空
  if (is.null(x = reduction)) {
  	# 则使用默认 Assay: "RNA"
    assay <- assay %||% DefaultAssay(object = object)
    data.use <- object[[assay]]
    cellnames.use <-  colnames(x = data.use) #细胞id
  } else {
  	# 否则，使用 reduction 对应的数据 
  	# > names(pbmc@reductions) #[1] "pca"  "umap" "tsne"
    data.use <- object[[reduction]]
    cellnames.use <- rownames(data.use)
  }
# > pbmc[["umap"]]
# A dimensional reduction object with key UMAP_ 
#  Number of dimensions: 2 
#  Projected dimensional reduction calculated:  FALSE 
#  Jackstraw run: FALSE 
#  Computed using assay: RNA 
# > pbmc[["RNA"]]
# Assay data with 13714 features for 2638 cells
# Top 10 variable features:
#  PPBP, LYZ, S100A9, IGLL5, GNLY, FTL, PF4, FTH1, GNG11, S100A8 


  #(A4) 根据分类名，获取旗下的 cell id。返回的是list()
  cells <- IdentsToCells(
    object = object,
    ident.1 = ident.1,
    ident.2 = ident.2,
    cellnames.use = cellnames.use
  )

  #(A5) 把倍数变化抛给了 FoldChange.Assay
  fc.results <- FoldChange(
    object = data.use,
    cells.1 = cells$cells.1,
    cells.2 = cells$cells.2,
    features = features,
    slot = slot,
    pseudocount.use = pseudocount.use,
    mean.fxn = mean.fxn,
    base = base,
    fc.name = fc.name
  )


  return(fc.results)
}










() FoldChange.Assay()
./seurat-4.1.0/R/differential_expression.R:1059:FoldChange.Assay <- function(


#' @importFrom Matrix rowMeans
#' @rdname FoldChange
#' @concept differential_expression
#' @export
#' @method FoldChange Assay
FoldChange.Assay <- function(
  object,
  cells.1,
  cells.2,
  features = NULL,
  slot = "data",
  pseudocount.use = 1,
  fc.name = NULL,
  mean.fxn = NULL,
  base = 2,
  ...
) {
  #(A1) 获取数据，等价于 pbmc@assays$RNA@data
  data <- GetAssayData(object = object, slot = slot)

  #(A2) 计算倍数变化使用的函数
  mean.fxn <- mean.fxn %||% switch(
    EXPR = slot,
    'data' = function(x) {
      return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
    },
    'scale.data' = rowMeans,
    function(x) {  # 这个没有名字的是默认分支
      return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
    }
  )

  #(A3) 如果底是e=2.71828，则不显示该值，否则显示
  # Omit the decimal value of e from the column name if base == exp(1)
  base.text <- ifelse(
    test = base == exp(1),
    yes = "",
    no = base
  )

  #(A4) 如果没有提供列名，则开始构建
  # 使用 scale.data 数据的，列名为 avg_diff
  # 否则，列名为 avg_log + 对数底 + FC，默认为 avg_log2FC
  fc.name <- fc.name %||% ifelse(
    test = slot == "scale.data",
    yes = "avg_diff",
    no = paste0("avg_log", base.text, "FC")
  )

  #(A5) 把计算任务推给 FoldChange.default
  FoldChange(
    object = data,
    cells.1 = cells.1,
    cells.2 = cells.2,
    features = features,
    mean.fxn = mean.fxn,
    fc.name = fc.name
  )
}






()FoldChange.default()
./seurat-4.1.0/R/differential_expression.R:1022:FoldChange.default <- function(

#' @param cells.1 Vector of cell names belonging to group 1
#' @param cells.2 Vector of cell names belonging to group 2

#' @param features Features to calculate fold change for.
#' If NULL, use all features

#' @importFrom Matrix rowSums
#' @rdname FoldChange
#' @concept differential_expression
#' @export
#' @method FoldChange default
FoldChange.default <- function(
  object,
  cells.1,
  cells.2,
  mean.fxn,
  fc.name,
  features = NULL,
  ...
) {
  #(A1) features 默认使用全部基因名
  features <- features %||% rownames(x = object)
  
  #(A2) 计算表达的百分比: Calculate percent expressed
  thresh.min <- 0

  #(A3) 只要大于0，就认为表达。
  # 统计一行中大于0的细胞个数，除以细胞总个数
  pct.1 <- round(
    x = rowSums(x = object[features, cells.1, drop = FALSE] > thresh.min) /
      length(x = cells.1),
    digits = 3
  )

  #(A4) 第2类的细胞表达百分比
  pct.2 <- round(
    x = rowSums(x = object[features, cells.2, drop = FALSE] > thresh.min) /
      length(x = cells.2),
    digits = 3
  )

  #(A5) 计算倍数变化 Calculate fold change
  # 使用对数公式: log(x/y)=log(x) - log(y)
  # 所以 log2( (m1+1)/(m2+1)) = log2(m1+1) - log2(m2+1)
  data.1 <- mean.fxn(object[features, cells.1, drop = FALSE])
  data.2 <- mean.fxn(object[features, cells.2, drop = FALSE])
  fc <- (data.1 - data.2)

  #(A6) 合并三列，并修改名字
  fc.results <- as.data.frame(x = cbind(fc, pct.1, pct.2))
  colnames(fc.results) <- c(fc.name, "pct.1", "pct.2")

  return(fc.results)
}










() FindMarkers.default()
../seurat-4.1.0/R/differential_expression.R:492:FindMarkers.default <- function(


#' @param cells.1 Vector of cell names belonging to group 1
#' @param cells.2 Vector of cell names belonging to group 2

#' @param counts Count matrix if using scale.data for DE tests. This is used for
#' computing pct.1 and pct.2 and for filtering features based on fraction
#' expressing
counts: 如果使用 scale.data 计算DEG，则使用counts计算 pct.1 和pct.2，及根据表达过滤基因

#' @param features Genes to test. Default is to use all genes

#' @param logfc.threshold Limit testing to genes which show, on average, at least
#' X-fold difference (log-scale) between the two groups of cells. Default is 0.25
#' Increasing logfc.threshold speeds up the function, but can miss weaker signals.
logfc.threshold: 限制做检验的基因数。默认是 0.25，越大越快，但会过滤掉越多的弱信号


#' @param test.use Denotes which test to use. Available options are: | 统计检验方法，默认是 wilcox
#' \itemize{
#'  \item{"wilcox"} : Identifies differentially expressed genes between two
#'  groups of cells using a Wilcoxon Rank Sum test (default)
#'  \item{"bimod"} : Likelihood-ratio test for single cell gene expression,
#'  (McDavid et al., Bioinformatics, 2013)
#'  \item{"roc"} : Identifies 'markers' of gene expression using ROC analysis.
#'  For each gene, evaluates (using AUC) a classifier built on that gene alone,
#'  to classify between two groups of cells. An AUC value of 1 means that
#'  expression values for this gene alone can perfectly classify the two
#'  groupings (i.e. Each of the cells in cells.1 exhibit a higher level than
#'  each of the cells in cells.2). An AUC value of 0 also means there is perfect
#'  classification, but in the other direction. A value of 0.5 implies that
#'  the gene has no predictive power to classify the two groups. Returns a
#'  'predictive power' (abs(AUC-0.5) * 2) ranked matrix of putative differentially
#'  expressed genes.
#'  \item{"t"} : Identify differentially expressed genes between two groups of
#'  cells using the Student's t-test.
#'  \item{"negbinom"} : Identifies differentially expressed genes between two
#'   groups of cells using a negative binomial generalized linear model.
#'   Use only for UMI-based datasets
#'  \item{"poisson"} : Identifies differentially expressed genes between two
#'   groups of cells using a poisson generalized linear model.
#'   Use only for UMI-based datasets
#'  \item{"LR"} : Uses a logistic regression framework to determine differentially
#'  expressed genes. Constructs a logistic regression model predicting group
#'  membership based on each feature individually and compares this to a null
#'  model with a likelihood ratio test.
#'  \item{"MAST"} : Identifies differentially expressed genes between two groups
#'  of cells using a hurdle model tailored to scRNA-seq data. Utilizes the MAST
#'  package to run the DE testing.
#'  \item{"DESeq2"} : Identifies differentially expressed genes between two groups
#'  of cells based on a model using DESeq2 which uses a negative binomial
#'  distribution (Love et al, Genome Biology, 2014).This test does not support
#'  pre-filtering of genes based on average difference (or percent detection rate)
#'  between cell groups. However, genes may be pre-filtered based on their
#'  minimum detection rate (min.pct) across both cell groups. To use this method,
#'  please install DESeq2, using the instructions at
#'  https://bioconductor.org/packages/release/bioc/html/DESeq2.html
#' }



#' @param min.pct  only test genes that are detected in a minimum fraction of
#' min.pct cells in either of the two populations. Meant to speed up the function
#' by not testing genes that are very infrequently expressed. Default is 0.1
min.pct: 在2群中至少一个有这个比例的表达量，才用于统计检验。
用于加速。主动跳过在2组中都不怎么表达的基因。默认是 0.1


#' @param min.diff.pct  only test genes that show a minimum difference in the
#' fraction of detection between the two groups. Set to -Inf by default
min.diff.pct: 2群中表达百分比的差异，至少要是该值，默认是 -Inf 


#' @param only.pos Only return positive markers (FALSE by default) 
only.pos: 是否仅返回正marker? 默认FALSE

#' @param verbose Print a progress bar once expression testing begins
verbose: 打印进度条

#' @param max.cells.per.ident Down sample each identity class to a max number.
#' Default is no downsampling. Not activated by default (set to Inf)
max.cells.per.ident: 每个类最多抽多少个样本。
默认是不抽样。

#' @param random.seed Random seed for downsampling
random.seed: 抽样的种子


#' @param latent.vars Variables to test, used only when \code{test.use} is one of
#' 'LR', 'negbinom', 'poisson', or 'MAST'
latent.vars: 统计检验的变量，只在 test.use 为 'LR', 'negbinom', 'poisson', or 'MAST' 时有用。
//todo 补充这几个统计的原理和方法。


#' @param min.cells.feature Minimum number of cells expressing the feature in at least one
#' of the two groups, currently only used for poisson and negative binomial tests
min.cells.feature: 2组中至少一组表达该基因的细胞最小个数。目前仅用于 泊松检验和 负二项分布检验
//todo 补充原理

#' @param min.cells.group Minimum number of cells in one of the groups
min.cells.group: 每组中最小细胞数。


#' @param pseudocount.use Pseudocount to add to averaged expression values when
#' calculating logFC. 1 by default.
pseudocount.use: 求log前加一个数，默认是1

#' @param fc.results data.frame from FoldChange
fc.results: 就是 函数 FoldChange() 的返回值，是一个3列数据框。

#' @param densify Convert the sparse matrix to a dense form before running the DE test. 
# This can provide speedups but might require higher memory; default is FALSE
densify: 做DE前，稀疏矩阵转为稠密矩阵。默认FALSE
这可以提速，但是占用更多内存。

#'
#'
#' @importFrom Matrix rowMeans
#' @importFrom stats p.adjust
#'
#' @rdname FindMarkers
#' @concept differential_expression
#' @export
#' @method FindMarkers default
#'
FindMarkers.default <- function(
  object,
  slot = "data",
  counts = numeric(),
  cells.1 = NULL,
  cells.2 = NULL,
  features = NULL,
  logfc.threshold = 0.25,
  test.use = 'wilcox',
  min.pct = 0.1,
  min.diff.pct = -Inf,
  verbose = TRUE,
  only.pos = FALSE,
  max.cells.per.ident = Inf,
  random.seed = 1,
  latent.vars = NULL,
  min.cells.feature = 3,
  min.cells.group = 3,
  pseudocount.use = 1,
  fc.results = NULL,
  densify = FALSE,
  ...
) {
  #(A1) 新函数: 见下文。
  ValidateCellGroups(
    object = object,
    cells.1 = cells.1,
    cells.2 = cells.2,
    min.cells.group = min.cells.group
  )

  #(A2) features 默认使用全部基因名
  features <- features %||% rownames(x = object)

  #(A3) 对于不需要过滤基因的DEG检测方法: > Seurat:::DEmethods_noprefilter() #[1] "DESeq2"
  # reset parameters so no feature filtering is performed
  if (test.use %in% DEmethods_noprefilter()) {
    features <- rownames(x = object) #使用全部基因
    min.diff.pct <- -Inf #组间差异百分比：-Inf
    logfc.threshold <- 0 #logfc 阈值: 0
  }

  #(A4) 数据使用哪个 slot?
  # 默认是 data, 就是 传入的矩阵obj
  # 如果是'scale.data'，则使用 counts: 做过滤
  data <- switch(
    EXPR = slot,
    'scale.data' = counts,
    object
  )

  #(A5) 选择基因: 基于表达百分比 feature selection (based on percentages)
  # 在两个群中，每个基因较大的表达百分比
  alpha.min <- pmax(fc.results$pct.1, fc.results$pct.2)
  # 命名向量：基因为名字
  names(x = alpha.min) <- rownames(x = fc.results)
  # 过滤掉在两者中表达百分比都很低的基因，返回至少在一个中pct超过阈值的基因名
  features <- names(x = which(x = alpha.min >= min.pct))
  # 如果基因长度为0，则警告：返回空df
  if (length(x = features) == 0) {
    warning("No features pass min.pct threshold; returning empty data.frame")
    return(fc.results[features, ])
  }

  #(A6) 最大百分比 - 最小百分比 = 基因在2组中表达百分比的差
  alpha.diff <- alpha.min - pmin(fc.results$pct.1, fc.results$pct.2)
  # 保留大于 最小差异 的基因名字
  features <- names(
    x = which(x = alpha.min >= min.pct & alpha.diff >= min.diff.pct)
  )
  # 如果基因长度为0，则警告：返回空df
  if (length(x = features) == 0) {
    warning("No features pass min.diff.pct threshold; returning empty data.frame")
    return(fc.results[features, ])
  }

  #(A7) 选择基因: 基于 logFC # feature selection (based on logFC)
  # 默认 slot="data"
  if (slot != "scale.data") {
  	# 获取第1列，就是 log2FC
    total.diff <- fc.results[, 1] #first column is logFC
    # 给vector加上基因名name
    names(total.diff) <- rownames(fc.results)
    # 如果只保留 正marker
    features.diff <- if (only.pos) {
      # 只保留 log2FC >=阈值 的基因名
      names(x = which(x = total.diff >= logfc.threshold))
    } else {
      # 否则，保留 |log2FC| >=阈值 的基因名
      names(x = which(x = abs(x = total.diff) >= logfc.threshold))
    }
    # 求前面 pct 过滤后的基因，和 logFC 过滤后的基因的交集
    features <- intersect(x = features, y = features.diff)
    # 如果长度为0，则警告：返回空df
    if (length(x = features) == 0) {
      warning("No features pass logfc.threshold threshold; returning empty data.frame")
      return(fc.results[features, ])
    }
  }

  #(A8) 如果细胞太多，则抽样取子集
  # subsample cell groups if they are too large
  # 如果设置过 max.cells.per.ident，则要抽样
  if (max.cells.per.ident < Inf) {
  	# 设置随机数种子
    set.seed(seed = random.seed)
    # 抽样
    if (length(x = cells.1) > max.cells.per.ident) {
      cells.1 <- sample(x = cells.1, size = max.cells.per.ident)
    }
    if (length(x = cells.2) > max.cells.per.ident) {
      cells.2 <- sample(x = cells.2, size = max.cells.per.ident)
    }
    # 如果 latent.vars 非空
    if (!is.null(x = latent.vars)) {
      # 则取出这几个细胞对应的子集
      latent.vars <- latent.vars[c(cells.1, cells.2), , drop = FALSE]
    }
  }

  #(A9) 做差异分析 PerformDE，见下文
  # 返回一个数据框：单列 p-value, 行名是基因名
  de.results <- PerformDE(
    object = object,
    cells.1 = cells.1,
    cells.2 = cells.2,
    features = features,
    test.use = test.use,
    verbose = verbose,
    min.cells.feature = min.cells.feature,
    latent.vars = latent.vars,
    densify = densify,
    ...
  )

  #(A10) 按列合并成df: p值，log2fc, pct.1, pct.2
  de.results <- cbind(de.results, fc.results[rownames(x = de.results), , drop = FALSE])
  
  #(A11) 如果只保留正marker：按照 log2FC >0 过滤
  if (only.pos) {
    de.results <- de.results[de.results[, 2] > 0, , drop = FALSE]
  }

  #(A12) 不用p值矫正的算法 "roc"
  if (test.use %in% DEmethods_nocorrect()) {
    de.results <- de.results[order(-de.results$power, -de.results[, 1]), ]
  } else {
  	# 需要做p值矫正的算法
  	# 按p-value排序，按 log2FC 倒序排列
    de.results <- de.results[order(de.results$p_val, -de.results[, 1]), ]
    # 做p值矫正: 使用 "bonferroni" 法，一个严格的方法。
    # 因为作者过滤掉的很多低表达基因没有做检验，它们p值按1算
    de.results$p_val_adj = p.adjust(
      p = de.results$p_val,
      method = "bonferroni",
      n = nrow(x = object)
    )
  }

  return(de.results)
}











() ValidateCellGroups()
$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -n "^\s*ValidateCellGroups" 2>/dev/null --color=auto
./seurat-4.1.0/R/differential_expression.R:2195:ValidateCellGroups <- function(

# FindMarkers helper function for cell grouping error checking
ValidateCellGroups <- function(
  object,
  cells.1,
  cells.2,
  min.cells.group
) {
  #(A1) 如果 cell id vector 长度为0，报错
  if (length(x = cells.1) == 0) {
    stop("Cell group 1 is empty - no cells with identity class ", cells.1)
  } else if (length(x = cells.2) == 0) {
    stop("Cell group 2 is empty - no cells with identity class ", cells.2)
    return(NULL) #这一句永远不会执行到：stop()程序已经退出
  #(A2) 如果有一组细胞数小于3个，则报错
  } else if (length(x = cells.1) < min.cells.group) {
    stop("Cell group 1 has fewer than ", min.cells.group, " cells")
  } else if (length(x = cells.2) < min.cells.group) {
    stop("Cell group 2 has fewer than ", min.cells.group, " cells")
  #(A3) 如果有 cell id 不在 obj 中，则报错
  } else if (any(!cells.1 %in% colnames(x = object))) {
    bad.cells <- colnames(x = object)[which(x = !as.character(x = cells.1) %in% colnames(x = object))]
    stop(
      "The following cell names provided to cells.1 are not present: ",
      paste(bad.cells, collapse = ", ")
    )
  } else if (any(!cells.2 %in% colnames(x = object))) {
    bad.cells <- colnames(x = object)[which(x = !as.character(x = cells.2) %in% colnames(x = object))]
    stop(
      "The following cell names provided to cells.2 are not present: ",
      paste(bad.cells, collapse = ", ")
    )
  }
}











()PerformDE()
$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -n "^\s*PerformDE" 2>/dev/null --color=auto
./seurat-4.1.0/R/differential_expression.R:1885:PerformDE <- function(

求P值。


PerformDE <- function(
  object,
  cells.1,
  cells.2,
  features,
  test.use,
  verbose,
  min.cells.feature,
  latent.vars,
  densify,
  ...
) {
  #(A1) 如果检验方法不需要 latent 变量，且提供了，则警告:
  if (!(test.use %in% DEmethods_latent()) && !is.null(x = latent.vars)) {
    warning(
      "'latent.vars' is only used for the following tests: ",
      paste(DEmethods_latent(), collapse=", "),
      call. = FALSE,
      immediate. = TRUE
    )
  }

  #(A2) 如果不是这3个，则需要检查...参数
  # > Seurat:::DEmethods_checkdots()
  #[1] "wilcox" "MAST"   "DESeq2"
  if (!test.use %in% DEmethods_checkdots()) {
    CheckDots(...)
  }

  #(A3) 数据框取子集
  data.use <- object[features, c(cells.1, cells.2), drop = FALSE]

  #(A4) 是否转为稠密矩阵
  if (densify){
    data.use <- as.matrix(x = data.use)
  }

  #(A5) 根据统计检验方法的名字，调用不同的函数
  de.results <- switch(
    EXPR = test.use,
    'wilcox' = WilcoxDETest( #默认是这个
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      verbose = verbose,
      ...
    ),

    'bimod' = DiffExpTest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      verbose = verbose
    ),

    'roc' = MarkerTest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      verbose = verbose
    ),

    't' = DiffTTest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      verbose = verbose
    ),

    'negbinom' = GLMDETest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      min.cells = min.cells.feature,
      latent.vars = latent.vars,
      test.use = test.use,
      verbose = verbose
    ),

    'poisson' = GLMDETest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      min.cells = min.cells.feature,
      latent.vars = latent.vars,
      test.use = test.use,
      verbose = verbose
    ),

    'MAST' = MASTDETest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      latent.vars = latent.vars,
      verbose = verbose,
      ...
    ),

    "DESeq2" = DESeq2DETest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      verbose = verbose,
      ...
    ),

    "LR" = LRDETest(
      data.use = data.use,
      cells.1 = cells.1,
      cells.2 = cells.2,
      latent.vars = latent.vars,
      verbose = verbose
    ),

    stop("Unknown test: ", test.use) #默认：不知道是哪个检验，报错
  )

  return(de.results)
}









() WilcoxDETest()
$ find . | grep "R$" | grep -v "testthat" | grep -v "notes" | xargs grep -n "^\s*WilcoxDETest" 2>/dev/null --color=auto
./seurat-4.1.0/R/differential_expression.R:2254:WilcoxDETest <- function(


秩和检验。使用 limma::rankSumTestWithCorrelation，作者认为这个实现更快。


# Differential expression using Wilcoxon Rank Sum
#
# Identifies differentially expressed genes between two groups of cells using
# a Wilcoxon Rank Sum test. Makes use of limma::rankSumTestWithCorrelation for a
# more efficient implementation of the wilcoxon test. Thanks to Yunshun Chen and
# Gordon Smyth for suggesting the limma implementation.
#
# @param data.use Data matrix to test
# @param cells.1 Group 1 cells
# @param cells.2 Group 2 cells
# @param verbose Print a progress bar
verbose: 是否打印进度条？

# @param ... Extra parameters passed to wilcox.test
#
# @return Returns a p-value ranked matrix of putative differentially expressed features
返回一个p-value 排序的可能差异表达的基因列表

#
#' @importFrom pbapply pbsapply  #带 pb=progressive bar 的 apply 函数
#' @importFrom stats wilcox.test
#' @importFrom future.apply future_sapply
#' @importFrom future nbrOfWorkers
#
# @export
#
# @examples
# data("pbmc_small")
# pbmc_small
# WilcoxDETest(pbmc_small, cells.1 = WhichCells(object = pbmc_small, idents = 1),
#             cells.2 = WhichCells(object = pbmc_small, idents = 2))
#
WilcoxDETest <- function(
  data.use,
  cells.1,
  cells.2,
  verbose = TRUE,
  ...
) {
  #(A1) 对 data 按细胞取子集
  data.use <- data.use[, c(cells.1, cells.2), drop = FALSE]
  #(A2) 生成和 cells.1 等长的向量，类似 1:n
  j <- seq_len(length.out = length(x = cells.1))
  
  #(A2) 决定使用哪个 apply 函数
  #见解析补充 - future 包
  my.sapply <- ifelse(
    #如果打印详情，且只有1个核,使用 pbsapply; 否则使用 future_sapply
    test = verbose && nbrOfWorkers() == 1, 
    yes = pbsapply, 
    no = future_sapply
  )

  #(A3) 检查是否溢出?? 
  overflow.check <- ifelse(
  	# 第一行 *  第一行，是na吗？ 是na则表示溢出内存了；不是na表示可以继续运算
    test = is.na(x = suppressWarnings(length(x = data.use[1, ]) * length(x = data.use[1, ]))),
    yes = FALSE,
    no = TRUE
  )

  #(A4) 检查R包是否正常安装: 见本文 3(8)，返回是一个 named vector
  limma.check <- PackageCheck("limma", error = FALSE)
  
  #(A5) 如果limma包安装了(T)，且没有溢出(T)
  if (limma.check[1] && overflow.check) {
  	#(B1) 则使用 limma::rankSumTestWithCorrelation 计算p值
    p_val <- my.sapply(   #默认是 pbsapply，它比 sapply 多了一个进度条
      X = 1:nrow(x = data.use),  #对每个基因编号循环
      FUN = function(x) {    #x是一个行号，一行一个基因
      	# 这个分组信号怎么传递的？
      	# 这个j从哪里来的？本函数开头第二行就是：j=1:length(cell.1)
        #return(min(2 * min(limma::rankSumTestWithCorrelation(index = j, statistics = data.use[x, ])), 1))
        return(
        	min(
        		2 * min(limma::rankSumTestWithCorrelation(index = j, statistics = data.use[x, ])), 
        		1)
        )
      }
    )
  } else {
    #(B2) 否则，使用R自带的 wilcox 计算p值
    # 先给出警告：建议安装 limma，及安装语句
    if (getOption('Seurat.limma.wilcox.msg', TRUE) && overflow.check) {
      message(
        "For a more efficient implementation of the Wilcoxon Rank Sum Test,",
        "\n(default method for FindMarkers) please install the limma package",
        "\n--------------------------------------------",
        "\ninstall.packages('BiocManager')",
        "\nBiocManager::install('limma')",
        "\n--------------------------------------------",
        "\nAfter installation of limma, Seurat will automatically use the more ",
        "\nefficient implementation (no further action necessary).",
        "\nThis message will be shown once per session"
      )
      options(Seurat.limma.wilcox.msg = FALSE)
    }
    
    # (C1) 新建只有行名的df
    group.info <- data.frame(row.names = c(cells.1, cells.2))
    # (C2)添加一列，根据行名填充分组
    group.info[cells.1, "group"] <- "Group1"
    group.info[cells.2, "group"] <- "Group2"
    # 分组列改为因子类型
    group.info[, "group"] <- factor(x = group.info[, "group"])
    # (C3) 对data.use 的列(cell id)排序
    data.use <- data.use[, rownames(x = group.info), drop = FALSE]
    # (C4) 求p值
    p_val <- my.sapply(
      X = 1:nrow(x = data.use),
      FUN = function(x) {
        return(wilcox.test(data.use[x, ] ~ group.info[, "group"], ...)$p.value)
        #根据我的测试，这个更快:
        # return(wilcox.test(data.use[x, cells.1], data.use[x, cells.2], ...)$p.value) 
      }
    )
  }
  # 返回数据框: 行名是基因名，单列 p-value
  return(data.frame(p_val, row.names = rownames(x = data.use)))
}
